home *** CD-ROM | disk | FTP | other *** search
- ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; The data in this file contains enhancments. ;;;;;
- ;;; ;;;;;
- ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
- ;;; All rights reserved ;;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; (c) Copyright 1980 Massachusetts Institute of Technology ;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
- (in-package "MAXIMA")
- (macsyma-module nalgfa)
-
- (declare-top(special alg-num vlist *var *nosplitf *algvar
- *denom *num *ans algfac* $nalgfac ALPHA)
- (genprefix nalg))
-
- (load-macsyma-macros rzmac ratmac)
-
- (setq alg-num 0)
-
- (defun new-alg nil
- (newsym (concat '$alg (setq alg-num (f1+ alg-num)))))
-
-
- (defun psqfrp (p var)
- (zerop (pdegree (pgcd p (pderivative p var)) var)))
-
- (defun pgsubst (val var p) ;;generalized psubst substitutes any
- (cond ((pcoefp p) p) ;;expression for any var in p
- ((eq var (car p))
- (cond ((pzerop val)
- (pterm (cdr p) 0))
- ((do ((ld (cadr p) (car a))
- (a (cdddr p) (cddr a))
- (ans (caddr p)
- (pplus
- (ptimes ans
- (pexpt val
- (f- ld (car a))))
- (cadr a))))
- ((null a) (ptimes ans (pexpt val ld)))))))
- ((pointergp var (car p)) p)
- ((do ((a (cdddr p) (cddr a))
- (ans (ptimes (list (car p) (cadr p) 1)
- (pgsubst val var (caddr p)))
- (pplus ans
- (ptimes (list (car p) (car a) 1)
- (pgsubst val var (cadr a))))))
- ((null a) ans)))))
-
- (defun pvsubst (nvar ovar p)
- (cond ((or (pcoefp p) (pointergp ovar (car p))) p)
- ((eq ovar (car p))
- (cons nvar (cdr p)))
- (t (pgsubst (make-poly nvar) ovar p))))
-
- (defun ordervar (var l)
- (let ((mvar (lmainvar l)))
- (cond ((null mvar) l)
- ((null (pointergp mvar var)) (cons var l))
- ((let ((newvar (gensym)))
- (setq genvar (append genvar (list newvar)))
- (setplist newvar (symbol-plist var))
- (set newvar (f1+ (symbol-value mvar)))
- (cons newvar
- (mapcar #'(lambda (p) (pvsubst newvar var p))
- l)))))))
-
- (defun lmainvar (l) ;;main var of list of poly's
- (do ((l l (cdr l))
- (v))
- ((null l) v)
- (cond ((pcoefp (car l)))
- ((null v) (setq v (caar l)))
- ((pointergp (caar l) v)
- (setq v (caar l))))))
-
- (defun presult (p1 p2 var) ;;change call in algsys?
- (let ((genvar genvar))
- (setq var (ordervar var (list p1 p2))
- p1 (cadr var)
- p2 (caddr var)
- var (car var))
- (cond ((zerop (pdegree p1 var))
- (cond ((zerop (pdegree p2 var)) 1)
- ((pexpt p1 (cadr p2)))))
- ((zerop (pdegree p2 var))
- (pexpt p2 (cadr p1)))
- ((resultant p1 p2)))))
-
- (defun pcoefvec (p)
- (cond ((pcoefp p) (list p))
- ((do ((l)
- (i (cadr p) (f1- i))
- (p (cdr p)))
- ((signp l i) (nreverse l))
- (push (cond ((and p (f= (car p) i))
- (prog1 (cadr p) (setq p (cddr p))))
- ( 0 ))
- l)))))
-
- (defun algtrace1 (bvec tvec)
- (do ((i (f- (length tvec) (length bvec)) (f1- i)))
- ((f= i 0) (algtrace* bvec tvec))
- (setq bvec (cons 0 bvec))))
-
- (defun algtrace* (bvec tvec)
- (do ((b bvec (cdr b))
- (tr (car (last bvec))
- (pplus tr (car (last b)))))
- ((null (cdr b)) tr)
- (or (pzerop (car b))
- (do ((l (cdr b) (cdr l))
- (tv tvec (cdr tv)))
- ((null l))
- (rplaca l (pdifference (car l)
- (ptimes (car b) (car tv))))))))
-
- (defun algtrace (r p)
- (cond ((eq (caar r) (car p))
- (ratreduce (algtrace1 (pcoefvec (car r))
- (cdr (pcoefvec p)))
- (cdr r)))
- ((ratreduce (pctimes (cadr p) (car r))
- (cdr r)))))
-
- (DEFMFUN $algtrace (r p var)
- (let ((varlist (list var))
- (genvar nil))
- (rdis* (algtrace (rform r) (car (rform p))))))
-
- (defun good-form (l) ;;bad -> good
- (do ((l l (cddr l))
- (ans))
- ((null l) (nreverse ans))
- (push (cons (cadr l) (car l)) ans)))
-
- (defun bad-form (l) ;;good -> bad
- (mapcar #'(lambda (q) (list (cdr q) (car q))) l))
-
- (DEFMFUN $algfac n
- (cond ((f= n 3) (apply '$pfactoralg (listify n)))
- (t (let ((varlist))
- (cond ((f= n 2)
- (newvar (arg 2))
- (cond ((alike1 (arg 2) (car varlist))
- ($pfactoralg (arg 1) nil (arg 2)))
- (t ($pfactoralg (arg 1)
- (arg 2)
- (car (last varlist))))))
- (t (newvar (arg 1))
- (setq varlist
- (mapcan #'(lambda (q) (cond ((algpget q)
- (list q))
- (t nil)))
- varlist))
- (cond ((f= (length varlist) 1)
- ($pfactoralg (arg 1)
- nil
- (car varlist)))
- ((f> (length varlist) 1)
- (merror "too many algebraics"))
- (t (merror "no algebraics")))))))))
-
- (DEFMFUN $pfactoralg (f p alg)
- (let ((varlist (list alg))
- (genvar) (vlist) (tellratlist) ($ratfac)
- ($gcd '$alGebraic)
- ($algebraic) ($ratalgdenom t)
- (*denom 1) (*num 1) (*ans))
- (cond ((and (null p) (radfunp alg t)) (newvar (cadr alg)))
- (t (newvar p)))
- (newvar1 f)
- (cond ((null vlist) (merror "attempt to factor a constant")))
- (setq varlist (nconc varlist (sortgreat vlist)))
- (cond (p (setq p (cadr (ratrep* p)))
- (push (cons alg (mapcar (function pdis) (cdr p)))
- tellratlist))
- (t (setq p (algpget alg))
- (setq p (pdifference
- (pexpt (cadr (ratrep* alg)) (car p))
- (cadr p)))))
- (setq $algebraic t)
- (setq f (cadr (ratrep* f)))
- (setq f (pfactoralg1 f p 0))
- (cons '(mtimes)
- (cons (rdis (ratreduce *num *denom))
- (mapcar 'pdis f)))))
-
- (declare-top(special adn*)) ;also redefine fact5 to call nalgfac correctly
-
- (defun nalgfac (p mp)
- (let ((*num 1) (*denom 1) (*ans) (algfac*) ($nalgfac)
- ($algebraic t) ($gcd '$algebraic))
- (setq p (pfactoralg1 p mp 0))
- (setq adn* (times adn* *denom))
- (cond ((equal *num 1) p)
- (t (cons *num p)))))
-
- (setq *nosplitf t)
-
- (defun pfactoralg1 (f p ind)
- (cond ((pcoefp f) (setq *num (ptimeschk f *num)) *ans)
- ((f= (cadr f) 1) (setq f (pshift f (car p) ind))
- (push (algnormal f) *ans)
- (setq f (rquotient f (car *ans))
- *denom (ptimeschk (cdr f) *denom)
- *num (ptimeschk (car f) *num))
- *ans)
- ((equal (cdr f) (cdr p))
- (push (pdifference (make-poly (car f)) (make-poly (car p))) *ans)
- (setq f (rquotient f (car *ans))
- *denom (ptimeschk (cdr f) *denom))
- (pfactoralg1 (car f) p ind))
- ((zerop (pdegree f (car p)))
- (mapc #'(lambda (q) (pfactoralg1 (pshift q (car p) -1) p (f1+ ind)))
- (let (($algebraic nil)
- ($gcd (car *gcdl*)))
- (pfactor1 f)))
- *ans)
- (t (do ((l (let (($algebraic nil)
- ($gcd (car *gcdl*)))
- (pfactor (algnorm f p)))
- (cddr l))
- (polys)
- (temp)
- (alg (car p)))
- ((null l)
- (setq *num (ptimeschk f *num))
- (mapc #'(lambda (q) (pfactoralg1
- (pshift q alg -1) p (f1+ ind)))
- polys)
- *ans)
- (cond ((pcoefp (car l)) nil)
- (t (setq temp (cond ((null (cddr l)) f)
- (t (pgcd f (car l)))))
- (cond ((f= (cadr temp) 1)
- (setq temp (algnormal temp))
- (push (pshift temp alg ind) *ans))
- ((f= (cadr l) 1)
- (setq temp (algnormal temp))
- (push (pshift temp alg ind) *ans)
- (or *nosplitf
- (setq *nosplitf
- (list (car l) temp ind))))
- (t (push temp polys)))
- (setq f (rquotient f temp)
- *denom (ptimeschk (cdr f) *denom)
- f (car f)))) ))))
-
- (defun pshift (f alg c)
- (if (= c 0) f
- (pgsubst (pplus (make-poly (car f)) (pctimes c (make-poly alg)))
- (car f) f)))
-
-
-
- (DEFMFUN $splitfield (p var)
- (let ((varlist)
- (genvar)
- (genpairs)
- (*algvar)
- ($gcd '$ALGEBRAIC))
- (newsym var)
- (setq *algvar (caar (new-alg)))
- (setq p (psplit-field (cadr (ratf p))))
- (cons
- '(mlist)
- (cons (pdis* (car p))
- (mapcar 'rdis* (cdr p))))))
-
- (defun psplit-field (p) ;modresult?
- (let ((l (mapcar #'(lambda (q) (psplit-field1 (cdr q)))
- (good-form (pfactor p)))) ;don't normalize lcfs?
- ($algebraic t))
- (if (null (cdr l)) (car l)
- (do ((l l (cdr l))
- (prim) (zeroes) (temp))
- ((null l) (cons (or prim '|$splits in q|) zeroes))
- (cond ((eq (caar l) 'linear)
- (setq zeroes (cons (cdar l) zeroes)))
- ((null prim)
- (setq prim (caar l)
- zeroes (nconc (cdar l) zeroes)))
- ((setq temp
- (primelmt (cons (car p) (cdr prim))
- (cons (car p) (cdaar l))
- *algvar)
- zeroes
- (nconc
- (mapcar
- #'(lambda (q)
- (ratgsubst (cadddr temp) (caaar l) q))
- (cdar l))
- (mapcar
- #'(lambda (q)
- (ratgsubst (caddr temp) (car prim) q))
- zeroes))
- prim (car temp))))))))
-
- (defun plsolve (p)
- (ratreduce (pterm (cdr p) 0)
- (pminus (pterm (cdr p) 1))))
-
-
- (defun psplit-field1 (p)
- ;;returns minimal poly and list of roots
- ;;p must be square free
- (*bind* ((minpoly (cons *algvar (cdr p)))
- (zeroes) ($algebraic t)
- ($ratalgdenom t))
- (if (equal (cadr p) 1) (return (cons 'linear (plsolve p))))
- (do ((polys (list p) )
- (nminpoly)
- (*nosplitf nil nil)
- (alpha (cons (make-poly (car minpoly)) 1)))
- ((null polys)
- (cons minpoly zeroes))
- (push alpha zeroes)
- (putprop (car minpoly) (cdr minpoly) 'tellrat)
- (rplaca polys
- (car
- (rquotient (pctimes (cdr alpha) (car polys))
- (pdifference
- (pctimes (cdr alpha) (pget (caar polys)))
- (car alpha)))))
- (setq polys
- (mapcan
- #'(lambda (q)
- (cond ((equal (cadr q) 1) ;;linear factor
- (push (plsolve q) zeroes)
- nil) ;;flush linear factor
- ((list q))))
- (mapcan #'(lambda (q)
- (let ((*ans) (*num 1) (*denom 1))
- (nreverse (pfactoralg1 q minpoly 0))))
- polys)))
- (when *nosplitf
- (setq nminpoly (car *nosplitf)
- *nosplitf (cdr *nosplitf))
- (putprop *algvar (cdr nminpoly) 'tellrat)
- (let ((beta
- (plsolve (pgcd (cons (caar *nosplitf) (cdr minpoly))
- (exchangevar (car *nosplitf) *algvar)))))
- (setq alpha (ratplus (cons (make-poly *algvar) 1)
- (rattimes (cons (f- (cadr *nosplitf)) 1)
- beta t)))
- (setq zeroes
- (mapcar
- #'(lambda (q) (ratgsubst beta (car minpoly) q))
- zeroes))
- (setq polys
- (mapcar
- #'(lambda (q) (car (rgsubst beta (car minpoly) q)))
- polys))
- (setq minpoly
- (cons *algvar (cdr nminpoly))))))))
-
- (defun exchangevar (poly var)
- (let ((newvar (gensym))
- (ovar (car poly)))
- (set newvar (f1+ (eval ovar)))
- (pvsubst ovar newvar
- (pvsubst var ovar
- (pvsubst newvar var poly)))))
-
- (defun rgsubst (val var p) ;;generalized psubst substitutes any
- (cond ((pcoefp p)
- (cons p 1)) ;;expression for any var in p
- ((eq var (car p))
- (cond ((pzerop val)
- (cons (pterm (cdr p) 0) 1))
- ((do ((ld (cadr p) (car a))
- (a (cdddr p) (cddr a))
- (ans (cons (caddr p) 1)
- (ratplus
- (rattimes ans
- (ratexpt val
- (f- ld (car a)))
- t)
- (cons (cadr a) 1))))
- ((null a) (rattimes ans (ratexpt val ld) t))))))
- ((pointergp var (car p)) (cons p 1))
- (t (let ((newsym (gensym)))
- (set newsym (f1+ (symbol-value (car p))))
- (rgsubst val newsym (pvsubst newsym var p))))))
-
- (defun ratgsubst (val var rat)
- (ratquotient (rgsubst val var (car rat))
- (rgsubst val var (cdr rat))))
-
- (defun algnorm (f p)
- (presult f p (car p)))
-
- (DEFMFUN $algnorm (r p var)
- (let ((varlist (list var))
- (genvar))
- (setq r (ratf r)
- p (cadr (ratf p)))
- (rdis* (cons (algnorm (cadr r) p)
- (algnorm (cddr r) p)))))
-
- (defun sqfrnorm (f p fvar) ;;f must be sqfr, p is minpoly, fvar # pvar
- (*bind* ((pvar (car p)))
- (setq f (cdr (ordervar pvar (list f p))) ;;new main var will be car of p
- p (cadr f) f (car f)) ;make mainvar of f = mainvar(p)
- (do ((i 0 (f1+ i))
- (dif (pdifference (make-poly fvar) (make-poly (car p))))
- (f f (pgsubst dif fvar f))
- (res))
- ((and (eq (car f) (car p))
- (setq res (primpart (algnorm f p)))
- (psqfrp res fvar))
- (list res
- (*bind* (($algebraic t)) ;;;modified f
- (putprop pvar (cdr p) 'tellrat)
- (pvsubst pvar (car p) f))
- (car p)
- p
- i)))))
-
- (defun primelmt (a b gvar &aux ($algebraic nil))
- ;;a is a poly with coeff's in k(b)
- ;;gvar is new variable
- (let ((norm (sqfrnorm (cons gvar (cdr a)) b gvar))
- (alpha) (beta) ($ratalgdenom t))
- (rplaca norm (primpart (car norm)))
- (putprop gvar (cdar norm) 'tellrat)
- (setq $algebraic t
- beta (pgcd (cadddr norm)
- (pvsubst (caddr norm)
- (car b)
- (cadr norm))))
- (setq beta (plsolve beta)
- alpha (ratplus (cons (make-poly gvar) 1)
- (rattimes (cons (f- (cadddr (cdr norm))) 1)
- beta t)))
- (list (car norm) ;;minimal poly
- (pplus (make-poly a) ;;new prim elm in old guys
- (list (car b) 1 (f- (cadddr (cdr norm)))))
- alpha beta))) ;;in terms of gamma
-
- (comment Discriminant of a Basis)
- (DEFMFUN $bdiscr n
- (let ((varlist) (genvar))
- (xcons (bdiscr (mapcar #'rform (listify (f1- n)))
- (car (rform (arg n))))
- (list 'mrat 'simp varlist genvar))))
-
- (defun bdiscr (l minp)
- (det
- (mapcar #'(lambda (q1)
- (mapcar #'(lambda (q2) (algtrace (ptimes q1 q2) minp)) l))
- l)))
-
- ;(declare-top(unspecial genvar varlist *var *num *ans alpha adn*))
-